      SUBROUTINE SCALEG (N,MA,A,MB,B,LOW,IGH,CSCALE,CPERM,WK)
CSE
C
CPS       PURPOSE:
C
C             SCALES THE MATRICES A AND B IN THE GENERALIZED
C             EIGENVALUE PROBLEM A*X = (LAMBDA)*B*X SUCH THAT THE
C             MAGNITUDES OF THE ELEMENTS OF THE SUBMATRICES OF A AND B
C             AS DETERMINED BY LOW AND IGH ARE AS CLOSE TO UNITY AS
C             POSSIBLE IN THE LEAST SQUARES SENSE.
C
CPE
CAS         ***  ARGUMENT LIST  ***
C
C  INPUT PARAMETERS:
C
C  N - THE ORDER OF THE MATRICES A AND B;
C
C  MA - THE ROW (FIRST) DIMENSION OF THE ARRAY A AS SPECIFIED
C       IN THE CALLING PROGRAM;
C
C  A - A REAL TWO-DIMENSIONAL ARRAY WITH ROW DIMENSION MA AND COLUMN
C      DIMENSION AT LEAST N CONTAINING THE MATRIX A AS THE PROBLEM
C      IS DEFINED IN THE FUNCTION STATEMENT ABOVE;
C
C  MB - THE ROW (FIRST) DIMENSION OF THE ARRAY B AS SPECIFIED IN
C       THE CALLING PROGRAM;
C
C  B - A REAL TWO-DIMENSIONAL ARRAY WITH ROW DIMENSION MB AND COLUMN
C      DIMENSION AT LEAST N CONTAINING THE MATRIX B AS THE PROBLEM
C      IS DEFINED IN THE FUNCTION STATEMENT ABOVE;
C
C  LOW - THE INTEGER SPECIFYING THE BEGINNING INDEX FOR THE ROWS AND
C        COLUMNS OF A AND B TO BE SCALED;
C
C  IGH - THE INTEGER SPECIFYING THE ENDING INDEX FOR THE ROWS AND
C        COLUMNS OF A AND B TO BE SCALED;
C
C  CPERM - AN ARRAY OF ANY DIMENSION THAT MUST BE PROVIDED BY THE
C          USER FOR WORK SPACE. CPERM MUST CONTAIN AT LEAST N
C          LOCATIONS. ONLY LOCATIONS LOW THROUGH IGH ARE REFERENCED;
C
C  WK - AN ARRAY OF ANY DIMENSION THAT MUST BE PROVIDED BY THE
C       USER FOR WORK SPACE. WK MUST CONTAIN AT LEAST 6*N LOCATIONS.
C       ONLY LOCATIONS LOW THROUGH IGH, N+LOW THROUGH N+IGH,...,
C       5*N+LOW THROUGH 5*N+IGH ARE REFERENCED.
C
C  OUTPUT PARAMETERS:
C
C  A - CONTAINS THE SCALED A MATRIX;
C
C  B - CONTAINS THE SCALED B MATRIX;
C
C  CSCALE - THE REAL ARRAY OF DIMENSION AT LEAST N CONTAINING IN
C           ITS LOW THROUGH IGH LOCATIONS THE INTEGER EXPONENTS OF 2
C           USED FOR THE COLUMN SCALING FACTORS. THE OTHER LOCATIONS
C           ARE NOT REFERENCED;
C
C  WK - CONTAINS IN ITS LOW THROUGH IGH LOCATIONS THE INTEGER
C       EXPONENTS OF 2 USED FOR THE ROW SCALING FACTORS.
C
CAE
C  REQUIRED SUBPROGRAMS - NONE
C
C  REQUIRED FUNCTIONS - ABS,SIGN,FLOAT,LOG10
C
C  REFERENCE - R.C. WARD, BALANCING THE GENERALIZED EIGENVALUE
C              PROBLEM, TECHNICAL REPORT ORNL/CSD-47, OAK RIDGE,
C              TENNESSEE, 1979.
C
C              [SAME TITLE] SIAM J. SCI. STAT. COMPUT., VOL 2, 
C                           NO 2, JUNE 1981, PP 141-152.
C
C ****
C
      DIMENSION A(MA,N),B(MB,N),CSCALE(N),CPERM(N),WK(N,6)
      DOUBLE PRECISION A,B,CSCALE,CPERM,WK,T,EW,EWC,FI,FJ,TA,
     1                 TB,TC,COR,SUM,BASL,BETA,CMAX,COEF,COEF2,
     2                 COEF5,ALPHA,GAMMA,PGAMMA
C
      IF (LOW .EQ. IGH) GO TO 405
      DO 210 I=LOW,IGH
         WK(I,1) = 0.D0
         CSCALE(I) = 0.D0
         WK(I,3) = 0.D0
         WK(I,4) = 0.D0
         WK(I,5) = 0.D0
         WK(I,6) = 0.D0
         WK(I,2) = 0.D0
         CPERM(I) = 0.D0
  210 CONTINUE
C ****
C     COMPUTE RIGHT SIDE VECTOR IN RESULTING LINEAR EQUATIONS
C ****
      BASL = DLOG10(2.D0)
      DO 240 I=LOW,IGH
         DO 240 J=LOW,IGH
            TB = B(I,J)
            TA = A(I,J)
            IF (TA .EQ. 0.D0) GO TO 220
            TA = DLOG10(DABS(TA)) / BASL
  220       IF (TB .EQ. 0.D0) GO TO 230
            TB = DLOG10(DABS(TB)) / BASL
  230       WK(I,5) = WK(I,5) - TA - TB
            WK(J,6) = WK(J,6) - TA - TB
  240 CONTINUE
C
      NR = IGH-LOW+1
      COEF = 1.D0 / FLOAT(2*NR)
      COEF2 = COEF * COEF
      COEF5 = .5D0 * COEF2
      NRP2 = NR + 2
      BETA = 0.D0
      IT = 1
C ****
C     START GENERALIZED CONJUGATE GRADIENT ITERATION
C ****
  250 CONTINUE
      EW = 0.D0
      EWC = 0.D0
      GAMMA = 0.D0
      DO 260 I=LOW,IGH
         GAMMA = GAMMA + WK(I,5)*WK(I,5) + WK(I,6)*WK(I,6)
         EW = EW + WK(I,5)
         EWC = EWC + WK(I,6)
  260 CONTINUE
      GAMMA = COEF*GAMMA - COEF2*(EW**2 + EWC**2)
     1        - COEF5 * (EW - EWC)**2
      IF (IT .NE. 1) BETA = GAMMA / PGAMMA
      T = COEF5 * (EWC - 3.*EW)
      TC = COEF5 * (EW - 3.*EWC)
      DO 270 I=LOW,IGH
         WK(I,2) = BETA*WK(I,2) + COEF*WK(I,5) + T
         CPERM(I) = BETA*CPERM(I) + COEF*WK(I,6) + TC
  270 CONTINUE
C ****
C     APPLY MATRIX TO VECTOR
C ****
      DO 300 I=LOW,IGH
         KOUNT = 0
         SUM = 0.D0
         DO 290 J=LOW,IGH
            IF (A(I,J) .EQ. 0.D0) GO TO 280
            KOUNT = KOUNT + 1
            SUM = SUM + CPERM(J)
  280       IF (B(I,J) .EQ. 0.D0) GO TO 290
            KOUNT = KOUNT + 1
            SUM = SUM + CPERM(J)
  290    CONTINUE
         WK(I,3) = FLOAT(KOUNT)*WK(I,2) + SUM
  300 CONTINUE
      DO 330 J=LOW,IGH
         KOUNT = 0
         SUM = 0.D0
         DO  320 I=LOW,IGH
            IF (A(I,J) .EQ. 0.D0) GO TO 310
            KOUNT = KOUNT + 1
            SUM = SUM + WK(I,2)
  310       IF (B(I,J) .EQ. 0.D0) GO TO 320
            KOUNT = KOUNT + 1
            SUM = SUM + WK(I,2)
  320    CONTINUE
         WK(J,4) = FLOAT(KOUNT)*CPERM(J) + SUM
  330 CONTINUE
C
      SUM = 0.D0
      DO 340 I=LOW,IGH
         SUM = SUM + WK(I,2)*WK(I,3) + CPERM(I)*WK(I,4)
  340 CONTINUE
C
C     IF GAMMA=0.0 AND SUM=0.0 GO TO END OF ITERATION
      IF(GAMMA.EQ.0.0 .AND. SUM.EQ.0.0) GO TO 370
C
      ALPHA = 0.0
      ALPHA = GAMMA / SUM
C ****
C     DETERMINE CORRECTION TO CURRENT ITERATE
C ****
      CMAX = 0.D0
      DO 350 I=LOW,IGH
         COR = ALPHA * WK(I,2)
         IF (DABS(COR) .GT. CMAX) CMAX = DABS(COR)
         WK(I,1) = WK(I,1) + COR
         COR = ALPHA * CPERM(I)
         IF (DABS(COR) .GT. CMAX) CMAX = DABS(COR)
         CSCALE(I) = CSCALE(I) + COR
  350 CONTINUE
      IF (CMAX .LT. .5D0) GO TO 370
      DO 360 I=LOW,IGH
         WK(I,5) = WK(I,5) - ALPHA * WK(I,3)
         WK(I,6) = WK(I,6) - ALPHA * WK(I,4)
  360 CONTINUE
      PGAMMA = GAMMA
      IT = IT + 1
      IF (IT .LE. NRP2) GO TO 250
C ****
C     END GENERALIZED CONJUGATE GRADIENT ITERATION
C ****
  370 DO 380 I=LOW,IGH
         IR = WK(I,1) + DSIGN(.5D0,WK(I,1))
         WK(I,1) = IR
         JC = CSCALE(I) + DSIGN(.5D0,CSCALE(I))
         CSCALE(I) = JC
  380 CONTINUE
C ****
C     SCALE A AND B
C ****
      DO 400 I=1,IGH
         IR = WK(I,1)
         FI = 2.D0**IR
         IF (I .LT. LOW) FI = 1.D0
         DO 400 J=LOW,N
            JC = CSCALE(J)
            FJ = 2.D0**JC
            IF (J .LE. IGH) GO TO 390
            IF (I .LT. LOW) GO TO 400
            FJ = 1.D0
  390       A(I,J) = A(I,J)*FI*FJ
            B(I,J) = B(I,J)*FI*FJ
  400 CONTINUE
  405 CONTINUE
C ****
C     EXIT SCALEG
C ****
      RETURN
      END
